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OBJECTIVES 


The  objectives  of  this  work  are  to  develop  and  implement  numerical  tech¬ 
niques  that  will  enable  us  to  perform  simulations  of  spatially  evolving 
phenomena.  In  particular,  we  are  interested  in  studying  the  phenomenon  of 


turbulent  diffusion  flame  lift-off. 


The  flame  lift-off  phenomenon  occurs  when  the  speed  of  the  fuel  jet  exceeds 
a  certain  value.  The  flame  detaches  from  the  exit  and  is  stabilized  at  a  cer¬ 
tain  distance  downstream.  Different  mechanisms  have  been  proposed  to  explain 
the  lift-off  phenomenon.  Generally,  it  is  believed  that  the  phenomenon  is 
similar  to  that  in  premixed  gas  combustion  (Brzustowski ,  1980;  Gunther  et  al., 
1981).  At  the  jet  exit,  the  velocity  of  the  fluid  in  the  mixing  layer  is 
higher  than  the  flame  propagation  speed.  The  flame  cannot  be  stabilized  at 
that  location.  As  the  mixing  proceeds,  the  velocity  of  the  fluid  in  the  mixing 
layer,  where  the  species  ratio  roughly  becomes  stoichiometric,  decreases  to  the 
flame  speed.  The  flame  is  stabilized  at  that  position.  This  model  assumes 
that  the  mixing  of  the  fuel  and  the  oxidizer  as  predicted  by  the  nonreacting 
turbulence  model  reaches  the  molecular  level  everywhere  in  the  mixing  layer. 
Direct  numerical  simulations  of  a  reacting  mixing  layer  (Riley  et  al.,  1986) 
indicated  that  this  may  not  be  the  case. 

Another  model  was  suggested  by  Peters  and  Williams  (1983).  They  proposed 
that,  in  a  turbulent  mixing  layer,  the  turbulent  eddies  produce  highly 
stretched  and  contorted  sheets  across  which  molecular  diffusion  of  species 
occurs.  Combustion  appears  as  laminar  flamelets  on  these  sheets.  Based  on 
the  theory  that  a  laminar  flame  will  extinguish  itself  if  the  strain  rate  is 
sufficiently  high  that  the  local  reduced  Damkohler  number  is  lower  than  a 
critical  value  (Linan,  1974),  Peters  and  Williams  predicted  that  the  strain 
rate  near  the  exit  of  the  jet  is  too  high  for  the  flame  to  exist.  As  the 
mixing  layer  grows,  the  strain  rate  is  reduced  to  a  sufficiently  low  level 
that  a  flame  can  be  sustained.  They  proceeded  to  estimate  the  strain  rate  in 
the  mixing  layer.  Through  a  statistical  theory  that  predicts  disruption  of  a 
continuous  sheet  due  to  "holes"  on  the  sheet,  they  were  able  to  obtain 
quantitative  results  for  lift-off  distance.  Their  theoretical  predictions  of 
lift-off  heights  are  of  the  right  "order  of  magnitude"  for  limited  methane 


flame  data. 
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The  results  of  Peters  and  Williams,  based  on  simple  turbulence  models,  are 
encouraging.  However,  a  detailed  numerical  simulation  of  the  interaction 
between  fluid  dynamical  and  combustion  processes  that  occur  in  diffusion 
flames  would  be  invaluable  to  further  our  understanding  of  the  physics. 

Direct  numerical  simulations  consist  of  accurately  solving  appropriate 
convection-diffusion-reaction  transport  equations  by  means  of  very  accurate 
numerical  algorithms  so  that  no  turbulence  modeling  is  required.  The  direct 
numerical  simulation  technique  has  recently  been  successfully  applied  to 
chemically  reacting  flows.  Riley  et  al.  (1986)  considered  the  three- 
dimensional  temporally  growing  mixing  layer  under  the  simplest  possible 
assumption  of  a  constant-rate  chemical  reaction  with  no  heat  release.  The 
main  contribution  of  this  work  is  the  understanding  of  the  effects  of  three- 
dimensional  mixing  and  diffusion  of  the  species  on  the  chemical  reaction. 
McMurtry  et  al.  (1985)  considered  the  effects  of  chemical  heat  release  and  the 
resulting  density  variation  on  the  fluid  motion  for  a  two-dimensional  mixing 
layer.  The  fluid  dynamics  and  the  chemical  reaction  are  truly  coupled  in  this 
work,  and  the  interplay  between  the  two  are  discussed.  However,  the  assumption 
of  a  constant  chemical  reaction  is  still  employed.  We  intend  to  extend  the 
simulations  to  study  the  local  quenching  and  the  lift-off  of  a  diffusion  flame. 
The  existing  methods  must  be  modified  to  include  (1)  a  temperature-dependent 
chemical  reaction  rate  and  (2)  a  spatially  developing  flow. 
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STATUS  OF  THE  RESEARCH 

In  the  first  ”^ar  of  this  research,  we  focused  on  an  initial  understanding 
of  the  physical  aspects  of  the  problem  through  simulations  of  two-dimensional 
flows.  During  the  second  year,  we  mainly  concentrated  on  construction  of  an 
accurate  three-dimensional  numerical  code. 

In  the  first  year,  we  examined  the  effects  of  large  coherent  structures  in 
two-dimensional,  unpremixed  chemically  reacting  mixing  layers  under  both 
temporally  evolving  and  spatially  developing  assumptions.  The  results  were 
reported  in  detail  by  Givi  et  al.  (1986)  and  are  only  summarized  here. 

(1)  In  a  two-dimensional  temporally  evolving  mixing  layer,  a  temperature- 
dependent  chemical  reaction  was  incorporated  into  a  computer  code  that 
uses  pseudospectral  numerical  methods.  The  nonequilibrium  effects  leading 
to  the  local  quenching  of  a  diffusion  flame  were  investigated.  The  results 
indicated  that  the  most  important  parameter  to  be  considered  for  flame 
extinction  is  the  local  instantaneous  scalar  dissipation  rate  conditioned 
at  the  scalar  stoichiometric  value.  At  locations  where  this  value  is 
increased  beyond  a  critical  value,  the  local  temperature  decreases  and  the 
instantaneous  reaction  rate  drops  to  zero,  leading  to  the  local  quenching 
of  the  flame.  This  work  was  published  by  Givi  et  al.  (1987). 

(2)  For  the  purpose  of  simulating  spatially  developing  flows,  a  two-dimensional 
hybrid  pseudospectral-f inite  difference  code  was  developed.  The  resulting 
code  was  used  for  the  simulation  of  the  pre-transitional  region  of  a 
laboratory  mixing  layer.  The  asynmetric  nature  of  the  mixing  process  was 
numerically  simulated  and  the  "preferred"  mixed  concentration  value 
(Masutani  and  Bowman,  1987)  was  numerically  calculated  by  constructing  the 
profiles  of  the  probability  density  function  of  a  passive  scalar  quantity 
across  the  shear  layer.  The  results  of  this  simulation  were  also  used  to 
explain  the  shortcomings  associated  with  turbulence  models  using  gradient 
diffusion  approximations  to  model  the  turbulent  flux  of  the  pdf,  such  as 
the  one  used  previously  by  Givi  et  al.  (1985).  This  work  was  published  by 
Givi  and  Jou  (1987). 

To  understand  the  exact  mechanism  of  the  lift-off  in  turbulent  flames, 
direct  numerical  simulations  of  three-dimensional,  spatially  evolving  flows 
under  the  influence  of  finite-rate  chemical  reactions  are  required.  For  this 
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purpose,  during  the  second  year  we  have  concentrated  on  constructing  an 
accurate  numerical  code  for  the  simulations  of  such  flows.  In  this  code,  the 
governing  equations  describing  the  hydrodynamic  variables  (velocities  and 
pressure)  and  scalar  quantities  (such  as  concentrations  and  mixture  tempera¬ 
ture)  are  solved  by  means  of  a  spectral  element  method.  This  method  combines 
the  versatility  of  finite  element  techniques  with  the  accuracy  of  spectral 
methods  in  a  more  flexible  manner  than  that  found  in  either  technique  alone 
(Patera,  1984).  This  approach  permits  arbitrary  boundary  conditions, 
therefore  allowing  us  to  simulate  spatially  evolving  flows  more  accurately. 

The  geometrical  configuration  of  the  type  of  the  flow  to  be  considered  is 
presented  in  Fig.  1.  The  flow  is  three-dimensional  and  evolves  spatially  in 
the  streamwise  (x)  direction.  An  explanation  of  the  spectral  element  algorithm 
is  given  in  Appendix  I,  and  a  detailed  explanation  of  the  computer  code  is 
Given  by  Givi  and  Israeli  (1987).  Here,  we  summarize  the  properties  of  the 
technique . 

(1)  The  streamwise  direction  is  divided  into  (NB)  number  of  blocks.  At  any 
time  during  computation  only  one  block  is  in  fast  memory  core  while  the 
other  blocks  are  in  secondary  storage  (such  as  the  SSD  or  disk)  and  can  be 
accessed  when  required.  Division  of  the  domain  in  this  manner  has  been 
done  to  ensure  portability  of  the  code  for  different  computers,  especially 
those  with  builted  fast  memory.  For  example,  on  a  CRAY-2  computer  (with 
256M  words  of  memory)  NB  can  be  set  equal  to  1  with  a  minimum  amount  of 
data  buffering,  and  the  values  of  NB  can  be  increased  to  accommodate 
machines  with  smaller  memories. 

(2)  Each  block  is  divided  into  NE  number  of  e lements ,  and  each  element  consists 
of  NI  number  of  Chebyshev  collocation  spectral  modes  in  the  streamwise 
direction.  This  division  has  the  advantage  that,  by  increasing  the  number 
of  elements  (increasing  NE)  for  better  resolution,  fewer  Chebyshev  modes 
(NI)  are  required.  Therefore,  the  distance  between  two  neighboring 
Chebyshev  points,  particularly  near  the  boundaries  of  the  elements,  will 
not  be  unreasonably  small.  Thus,  the  time  stepping  stability  constraints 
can  be  reduced . 

(3)  The  flow  is  assumed  periodic  in  the  spanwise  direction  (z),  and  impermeable 
free-slip  boundary  conditions  are  employed  in  the  cross-stream  direction 
(y).  Pseudospectral  methods  with  Fourier  transforms  are  used  in  these 
directions . 
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(4)  The  transport  equations  are  solved  in  Fourier  space,  after  the  Fourier 

transformation  of  the  variables  in  the  y  and  z  directions.  The  evaluation 
of  the  nonlinear  convection  terms,  however,  is  performed  in  physical 
space.  A  second-order-accurate  Adams -Bash for th  technique  (Roach,  1972)  is 
employed  for  the  time  discretization. 

The  hybrid  spectral  element  algorithm  used  in  the  streamwise  direction  has  the 
advantage  that  the  accuracy  of  spectral  methods  is  maintained  in  the  flow 
direction.  In  most  previous  work,  the  streamwise  derivatives  were  usually 
discretized  by  second-order-accurate  finite  difference  methods,  where  the 
convective  terms  were  approximated  by  an  upwind  difference  scheme  (Givi  and 
Jou,  1987;  Lowery  et  al.,  1987).  Second-order  central  differencing  requires  a 
maximum  cell  Reynolds  (Peclet)  number  of  2  (Roach,  1972),  which  can  require  an 
excessive  number  of  grid  points  in  the  x-direction.  The  upwind  scheme  results 
in  an  improvement  of  the  capability  in  resolving  sharp  gradients  for  moderate 
Reynolds  number  flow  simulations.  However,  it  results  in  the  addition  of  an 
artificial  numerical  viscosity  to  the  fluid  viscosity,  which  is  not  desirable, 
particularly  in  diffusion  flame  simulations.  The  implementation  of  the 
spectral  element  methods  in  the  streamwise  direction  will  allow  us  to  simulate 
flows  with  higher  Reynolds  numbers  accurately  without  the  addition  of  any 
artificial  viscosity. 

We  have  just  completed  the  construction  and  debugging  of  the  code.  The 
finite  element  matrices  and  the  Chebyshev  discretization  routines  are 
accurately  constructed.  Some  simple  problems  with  known  analytical  solutions 
have  been  used  as  test  problems.  The  results  of  preliminary  tests  using 
spectral  element  methods  show  excellent  agreement  with  the  analytical 
solutions  (see  Appendix  I). 

After  validating  the  code  by  comparison  of  the  data  with  known  analytical 
results,  simulations  of  a  laminar  three-dimensional  plane  jet  were  performed. 
In  Fig.  2  we  present  contour  plots  of  the  instantaneous  axial  velocity  of  the 
jet  in  a  x-y  plane.  The  figure  corresponds  to  a  time  when  the  fluid  has  swept 
the  computational  domain  once.  A  very  coarse  grid  was  chosen  (16**3)  for  the 
purpose  of  demonstration.  As  shown  on  the  figure,  the  linear  growth  of  the 
laminar  jet  is  well  predicted.  The  computational  time  for  this  simulation  was 
about  0.61  sec/time  step  on  a  CRAY-XMP,  which  is  very  encouraging  and 
indicates  that  higher  resolution  simulations  can  be  performed. 
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CURRENT  ACTIVITIES 

We  are  presently  starting  the  simulations  of  turbulent  diffusion  flames 
with  finite  rate  chemical  reactions.  The  resulting  code  will  be  used  to  study 
a  three-dimensional  turbulent  flow  field  in  a  spatially  evolving  jet  and  the 
structure  of  the  "holes"  on  the  flame  sheet.  Variable  density  effects  caused 
by  the  heat  release  will  then  be  added  to  the  code  at  a  later  stage. 

PUBLICATIONS 

The  following  written  manuscripts  have  resulted  from  our  efforts  during 
the  second  year: 

(1)  Givi,  P. ,  and  Jou,  W.-H. ,  "Mixing  and  Chemical  Reactions  in  a  Spatially 
Developing  Mixing  Layer,"  accepted  for  publication  in  Journal  of 
Non-Equilibrium  Thermodynamics,  to  appear  (1987)  ;  also  presented  at  Spring 
Technical  Meeting  of  The  Combustion  Institute,  Central  States  Section, 
Cleveland,  Ohio,  May  5-6,  1986. 

(2)  Givi,  P. ,  Jou,  W.-H.,  and  Metcalfe,  R.  W.,  "Flame  Extinction  in  a 
Temporally  Developing  Mixing  Layer,"  Proceedings  of  the  21st  Symposium 
(International)  on  Combustion,  to  appear  (1987). 

The  following  paper  was  in  part  a  result  of  this  work: 

(1)  McMurtry,  P.  A.,  and  Givi,  P. ,  ''Direct  Numerical  Simulations  of  Scalar 
Dissipation  in  a  Turbulent  Mixing  Layer",  SIAM  Conference  on  Numerical 
Combustion,  P.  A23,  Paper  #77,  San  Francisco,  CA,  March  9-11  (1987). 

INTERACTIONS 

In  addition  to  the  above-mentioned  papers,  parts  of  this  work  were 
presented  at  the  following  seminars: 

(1)  Givi,  P. ,  "Model  Free  Simulations  of  Turbulent  Diffusion  Flames,"  Invited 
Seminar,  Graduate  Seminar  Series,  Winter  Quarter,  Department  of  Mechanical 
Engineering,  University  of  Washington,  Seattle,  WA,  February  3,  1987. 

(2)  Givi,  P.,  "Direct  Numerical  Simulations  of  Turbulent  Diffusion  Jet  Flames," 
Invited  Seminar,  Fall  Seminar  Series,  Department  of  Mechanical  and 
Aerospace  Engineering,  Case-Western  Reserve  University,  Cleveland,  Ohio, 
September  26,  1986. 

(3)  Givi,  P. ,  "Fundamental  Studies  on  Turbulent  Combustion  and  Spray  Combustion, 
ICOMP,  NASA  Lewis  Research  Center,  Cleveland,  Ohio,  September  10,  1986. 
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APPENDIX  I 

A  SUMMARY  OF  THE  SPECTRAL  ELEMENT  METHOD 


FORMULATION 


The  geometrical  configuration  is  shown  in  Fig.  1.  The  flow  evolves 
spatially  in  the  streamwise  direction,  x,  is  periodic  in  the  spanwise 
direction,  z,  and  impermeable  boundary  conditions  are  used  in  the  cross- stream 
direction,  y.  For  the  purpose  of  demonstration,  only  the  normalized  transport 
equations  of  the  hydrodynamical  variables  (U  and  P)  are  considered: 


v-  U  - 


PU 


U  x-/2_  -  Vp  4.  J_  V  U 

Re 


where 


P  *  Dynamic  pressure  =  p  +  (1/2)U*U 
ft  =  Vorticity  *  V  x  U 
Re  -  Reynolds  mumber 


Employing  a  second-order  Adams-Bashfor th  scheme  for  the  temporal 
discretization,  we  would  have: 

.  *  ,  ,0  n  .  _  \H-t 


y-  U  -  3  /U  xSlf.  _L  fOx-^-) 
M:  -  -2 


W-  u') 


A*t 


J_  v*iT'-vP MJ 


where  At  is  the  computational  time  increment,  and  the  superscripts  n,  *,  and 
n+1  refer  to  the  previous,  intermediate,  and  future  times,  respectively. 
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Equation  (3)  constitutes  the  nonlinear  step,  and  the  derivatives  on  the 
RHS  of  this  equation  can  be  evaluated  by  the  use  of  Chebyshev  spectral  methods 
(Gottlieb  and  Orszag,  1977).  The  cross  terms  are  evaluated  in  physical  space, 
and  the  velocities  are  updated  in  time  by  the  second-order-accurate  Adams- 
Bashforth  scheme. 

Equations  (4)  and  (5)  constitute  the  pressure  step  and  the  viscous  step, 
respectively.  Their  numercial  simulations  can  be  achieved  by  finite  element 
methods . 

Taking  the  Fourier  transform  of  Equation  (5),  we  have: 

S  _  [(kV?)+||]  =  Re  VP*.  S|  U*  (6, 

In  this  equation,  -  indicates  the  Fourier  domain,  K  and  L  are  the  Fourier 
wave  numbers  in  the  y  and  z  directions,  respectively.  The  subscript  x  denotes 
the  derivative  in  the  streamwise  direction.  Following  the  same  procedure  for 
the  pressure  equation,  we  have:  / 


Pt  -  Ck  +  lx)  P 


y.  v 

bX- 


Note  that  Equations  (6)  and  (7)  can  be  written  in  a  generalized  form: 


Ux*  -  A  u  *■  ? 


2 

where  A  represents  the  coeffients  of  U  and  P,  and  f  represents  the  RHS. 

Employing  the  variational  principle  equivalent  to  the  solution  of  Equation 
(8),  we  have  the  following  elemental  equation: 


Ae  U  -  A*  y  =  Be| 


Where  Ae  and  Beware  (NI+1 )x(NI+l )  matrices,  and  NI  is  the  number  of  intervals 

in  an  element  and  also  the  degree  of  the  interpolating  polynominal.  For  a 

,2 

given  A  ,  we  have 


=  Ae  -  A  Be 


Ce  =  hS 
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Therefore, 


Ce  U  =  Be  £  =•  Pn 

Assume  U  and  ft  are  decomposed  as  shown  in  Fig.  3: 


where  the  subscripts  f,  I  and  L  indicate  the  first  element,  the 
elements,  and  the  last  element,  respectively. 

Also,  let  us  decompose  the  matrix  Ce  as  follows: 


(11) 

(12) 

(13) 

interior 

(14) 


Therefore,  we  have  the  elemental  equation: 

Pe  I  g?T  i  Be  I  U? 
fiF  A1  RU  U 


Be.  fy=- 


Solving  for  the  interior  elements,  U,  we  have: 

uE=  (aj)”'  L&  -  ^  u^-£s=  uL  J 


and  denoting  (AI)  by  A,  we  have: 


ur  =  A  r\  U*  f\  EE  -  UL  A  RL 


TR-40 1/03-86 


10 


Consider  the  first  internal  equation: 

Be  IU  RL  U,\  zDe  Uz+  RF  (J*  +  Be  U3  =  R  +  <1! 

But,  from  Equation  (17): 

U,1  =  A  (  ?,*-  U,  RF_  Uz  Pjr) 

Uz1-  A  {£*-  UzRF-U3RX) 

Substituing  Equations  (20)  and  (21)  in  (19),  we  have: 

(Be-RLA  gF)U,  +  (2-De.- RL 5  Rf  -  Rf  J  RJ- )  U*. 

+  (Be  -  Rf  A  RL)  U3  =-  RL  A  r!  +  -R_FJ\  Rz. 

Equation  (22)  can  be  generalized  by  substituting 


1- ^j-l 

2- *j 


(Be- rl  £ RF)  Uj_,  +  C2-^- Sk  5  EL)  r 

+i.ee.-  pf  a  EL)  ^j+>  =  +Ri->  Rf-£F  £Ej<2 

Note  that  Equation  (23)  forms  a  tridiagonal  system  of  equations  that  can  be 
solved  easily  when  the  values  at  the  interior  planes  are  known. 
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TEST  PROBLEM 


To  test  the  viscous  solver,  the  following  example  is  considered: 

0**  _  A*  U  “  ° 

(24) 

U(*-=o)  -  ' 

=  ° 

The  exact  solution  of  this  equation  is: 


Equation  (24)  was  solved  by  the  spectral  element  technique  with  three 
elements  (NE  **  3)  and  two  Chebyshev  interpolation  points  within  each  element 
(a  total  of  7  points  in  the  x  direction). 

In  Fig.  4,  the  results  of  the  numerical  solution  are  compared  with  the 
analytical  results.  As  can  be  seen,  the  agreement  is  excellent. 
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□  ft  -  0.  SPECTRAL  ELEMEN1 
O  fi  -  5.  SPECTRAL  ELEMENT 


